clear all
close all
clc

N=100;
D=zeros(N+1,2);
i=0;
for P=0:10000/N:10000
    
    i=i+1;
    D(i,1)=P;  
    D(i,2)=1-1.9722e-11*P^2.4279;
    
end

% syms p
% f=fittype('0.00001+(1-0.00001)*(1+(a*p)^n)^(-1+1/n)','independent',{'p'},'coefficients',{'a','n'});
% [cfun,rsquare]=fit(D(:,1),D(:,2),f,'Lower',[0,0.8985],'upper',[10000,1],'startPoint',[4000, 0.989])
% S=cfun(D(:,1));

fo = fitoptions('Method','NonlinearLeastSquares','Lower',[0,0.8985],'upper',[10000,1]);
ft = fittype('0.000001+(1-0.000001)*(1+(a*x)^n)^(-1+1/n)','problem','n','options',fo);
[curve2,gof2] = fit(D(:,1),D(:,2),ft,'problem',2.6)

S=curve2(D(:,1));
plot(D(:,1),D(:,2),'ro',D(:,1),S,'b-');

IrreducibleSaturation=0.0001;
SWCCa=5.253e-5;
MatricSuction=1e4;
SWCCn=2.6;
SWCCm=1-1/SWCCn ;
IrreducibleSaturation+(1.0D0-IrreducibleSaturation)*(1.0D0+(SWCCa*MatricSuction)^SWCCn)^(-SWCCm)